home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / b / b.lha / B / src / bint / b1nuA.c < prev    next >
Encoding:
C/C++ Source or Header  |  1988-11-24  |  6.4 KB  |  253 lines
  1. /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1985. */
  2.  
  3. /*
  4.   $Header: b1nuA.c,v 1.4 85/08/22 16:50:25 timo Exp $
  5. */
  6.  
  7. /* Approximate arithmetic */
  8.  
  9. #include "b.h"
  10. #include "b0con.h"
  11. #include "b1obj.h"
  12. #include "b1num.h"
  13. #include "b3err.h" /* For still_ok */
  14.  
  15. /*
  16. For various reasons, on some machines (notably the VAX), the range
  17. of the exponent is too small (ca. 1.7E38), and we cope with this by
  18. adding a second word which holds the exponent.
  19. However, on other machines (notably the IBM PC), the range is sufficient
  20. (ca. 1E300), and here we try to save as much code as possible by not
  21. doing our own exponent handling.  (To be fair, we also don't check
  22. certain error conditions, to save more code.)
  23. The difference is made by #defining EXT_RANGE (in b1num.h), meaning we
  24. have to EXTend the RANGE of the exponent.
  25. */
  26.  
  27. #ifdef EXT_RANGE
  28. Hidden struct real app_0_buf = {Num, 1, -1, 0.0, -(double)Maxint};
  29.     /* Exponent must be less than any realistic exponent! */
  30. #else !EXT_RANGE
  31. Hidden struct real app_0_buf = {Num, 1, -1, 0.0};
  32. #endif !EXT_RANGE
  33.  
  34. Visible real app_0 = &app_0_buf;
  35.  
  36. /*
  37.  * Build an approximate number.
  38.  */
  39.  
  40. Visible real mk_approx(frac, expo) double frac, expo; {
  41.     real u;
  42. #ifdef EXT_RANGE
  43.     expint e;
  44.     if (frac != 0) frac = frexp(frac, &e), expo += e;
  45.     if (frac == 0.5) frac = 1, --expo; /* Assert 0.5 < frac <= 1 */
  46.     if (frac == 0 || expo < -BIG) return (real) Copy(app_0);
  47.     if (expo > BIG) {
  48.         error(MESS(700, "approximate number too large"));
  49.         expo = BIG;
  50.     }
  51. #else !EXT_RANGE
  52.     if (frac == 0.0) return (real) Copy(app_0);
  53.     frac= ldexp(frac, (int)expo);
  54. #endif EXT_RANGE
  55.     u = (real) grab_num(-1);
  56.     Frac(u) = frac;
  57. #ifdef EXT_RANGE
  58.     Expo(u) = expo;
  59. #endif EXT_RANGE
  60.     return u;
  61. }
  62.  
  63. /*
  64.  * Approximate arithmetic.
  65.  */
  66.  
  67. Visible real app_sum(u, v) real u, v; {
  68. #ifdef EXT_RANGE
  69.     real w;
  70.     if (Expo(u) < Expo(v)) w = u, u = v, v = w;
  71.     if (Expo(v) - Expo(u) < Minexpo) return (real) Copy(u);
  72.     return mk_approx(Frac(u) + ldexp(Frac(v), (int)(Expo(v) - Expo(u))),
  73.         Expo(u));
  74. #else !EXT_RANGE
  75.     return mk_approx(Frac(u) + Frac(v), 0.0);
  76. #endif !EXT_RANGE
  77. }
  78.  
  79. Visible real app_diff(u, v) real u, v; {
  80. #ifdef EXT_RANGE
  81.     real w;
  82.     int sign = 1;
  83.     if (Expo(u) < Expo(v)) w = u, u = v, v = w, sign = -1;
  84.     if (Expo(v) - Expo(u) < Minexpo)
  85.         return sign < 0 ? app_neg(u) : (real) Copy(u);
  86.     return mk_approx(
  87.         sign * (Frac(u) - ldexp(Frac(v), (int)(Expo(v) - Expo(u)))),
  88.         Expo(u));
  89. #else !EXT_RANGE
  90.     return mk_approx(Frac(u) - Frac(v), 0.0);
  91. #endif !EXT_RANGE
  92. }
  93.  
  94. Visible real app_neg(u) real u; {
  95.     return mk_approx(-Frac(u), Expo(u));
  96. }
  97.  
  98. Visible real app_prod(u, v) real u, v; {
  99.     return mk_approx(Frac(u) * Frac(v), Expo(u) + Expo(v));
  100. }
  101.  
  102. Visible real app_quot(u, v) real u, v; {
  103.     if (Frac(v) == 0.0) {
  104.         error(MESS(701, "in u/v, v is zero"));
  105.         return (real) Copy(u);
  106.     }
  107.     return mk_approx(Frac(u) / Frac(v), Expo(u) - Expo(v));
  108. }
  109.  
  110. /*
  111.     YIELD log"(frac, expo):
  112.         CHECK frac > 0
  113.         RETURN normalize"(expo*logtwo + log(frac), 0)
  114. */
  115.  
  116. Visible real app_log(v) real v; {
  117.     double frac = Frac(v), expo = Expo(v);
  118.     if (frac <= 0) {
  119.         error(MESS(702, "in log x, x is <= 0"));
  120.         return (real) Copy(app_0);
  121.     }
  122.     return mk_approx(expo*logtwo + log(frac), 0.0);
  123. }
  124.  
  125. /*
  126.     YIELD exp"(frac, expo):
  127.         IF expo < minexpo: RETURN zero"
  128.         WHILE expo < 0: PUT frac/2, expo+1 IN frac, expo
  129.         PUT exp frac IN f
  130.         PUT normalize"(f, 0) IN f, e
  131.         WHILE expo > 0:
  132.             PUT (f, e) prod" (f, e) IN f, e
  133.             PUT expo-1 IN expo
  134.         RETURN f, e
  135. */
  136.  
  137. Visible real app_exp(v) real v; {
  138. #ifdef EXT_RANGE
  139.     expint ei;
  140.     double frac = Frac(v), expo = Expo(v), new_expo;
  141.     static double canexp;
  142.     if (!canexp) canexp = floor(log(log(Maxreal)) / logtwo);
  143.     if (expo <= canexp) {
  144.         if (expo < Minexpo) return mk_approx(1.0, 0.0);
  145.         frac = ldexp(frac, (int)expo);
  146.         expo = 0;
  147.     }
  148.     else if (expo >= Maxexpo) {
  149.         /* Definitely too big (the real boundary is much smaller
  150.            but here we are in danger of overflowing new_expo
  151.            in the loop below) */
  152.         return mk_approx(1.0, Maxreal); /* Force an error! */
  153.     }
  154.     else {
  155.         frac = ldexp(frac, (int)canexp);
  156.         expo -= canexp;
  157.     }
  158.     frac = exp(frac);
  159.     new_expo = 0;
  160.     while (expo > 0 && frac != 0) {
  161.         frac = frexp(frac, &ei);
  162.         new_expo += ei;
  163.         frac *= frac;
  164.         new_expo += new_expo;
  165.         --expo;
  166.     }
  167.     return mk_approx(frac, new_expo);
  168. #else !EXT_RANGE
  169.     return mk_approx(exp(Frac(v)), 0.0);
  170. #endif !EXT_RANGE
  171. }
  172.  
  173. /*
  174.     YIELD (frac, expo) power" v":
  175.         \   (f*2**e)**v =
  176.         \ = f**v * 2**(e*v) =
  177.         \ = f**v * 2**((e*v) mod 1) * 2**floor(e*v) .
  178.         PUT exp"(v" prod" normalize"(log frac, 0)) IN temp1" \ = f**v
  179.         PUT expo*numval(v") IN ev \ = e*v
  180.         PUT exp(logtwo * (ev - floor ev))) IN temp2 \ = 2**(ev mod 1)
  181.         PUT temp1" IN f, e
  182.         RETURN normalize"(f*temp2, e + floor ev)
  183. */
  184.  
  185. Visible real app_power(u, v) real u, v; {
  186.     double frac = Frac(u);
  187. #ifdef EXT_RANGE
  188.     real logfrac, vlogfrac, result;
  189.     double expo = Expo(u), rest;
  190. #endif !EXT_RANGE
  191.     if (frac <= 0) {
  192.         if (frac < 0) error(MESS(703, "in 0**v, v is negative"));
  193.         if (v == app_0) return mk_approx(1.0, 0.0); /* 0**0 = 1 */
  194.         return (real) Copy(app_0); /* 0**x = 0 */
  195.     }
  196. #ifdef EXT_RANGE
  197.     frac *= 2, expo -= 1; /* Renormalize to 1 < frac <= 2, so log frac > 0 */
  198.     logfrac = mk_approx(log(frac), 0.0);
  199.     vlogfrac = app_prod(v, logfrac);
  200.     result = app_exp(vlogfrac);
  201.     /* But what if result overflows but expo is very negative??? */
  202.     if (still_ok) {
  203.         expo *= numval((value)v);
  204.         rest = expo - floor(expo);
  205.         frac = Frac(result) * exp(logtwo*rest);
  206.         expo = Expo(result) + floor(expo);
  207.     }
  208.     release((value)logfrac), release((value)vlogfrac), release((value)result);
  209.     return mk_approx(frac, expo);
  210. #else !EXT_RANGE
  211.     return mk_approx(exp(log(frac) * Frac(v)), 0.0);
  212. #endif !EXT_RANGE
  213. }
  214.  
  215. Visible int app_comp(u, v) real u, v; {
  216.     double xu, xv;
  217. #ifdef EXT_RANGE
  218.     double eu, ev;
  219. #endif EXT_RANGE
  220.     if (u == v) return 0;
  221.     xu = Frac(u), xv = Frac(v);
  222. #ifdef EXT_RANGE
  223.     if (xu*xv > 0) {
  224.         eu = Expo(u), ev = Expo(v);
  225.         if (eu < ev) return xu < 0 ? 1 : -1;
  226.         if (eu > ev) return xu < 0 ? -1 : 1;
  227.     }
  228. #endif EXT_RANGE
  229.     if (xu < xv) return -1;
  230.     if (xu > xv) return 1;
  231.     return 0;
  232. }
  233.  
  234. Visible value app_floor(u) real u; {
  235. #ifdef EXT_RANGE
  236.     integer v, w;
  237.     value twotow, result;
  238.     if (Expo(u) <= Dblbits)
  239.         return (value)
  240.             mk_int(floor(ldexp(Frac(u),
  241.                 (int)(Expo(u) < 0 ? -1 : Expo(u)))));
  242.     v = mk_int(ldexp(Frac(u), Dblbits));
  243.     w = mk_int(Expo(u) - Dblbits);
  244.     twotow = power((value)int_2, (value)w);
  245.     result = prod((value)v, twotow);
  246.     release((value) v), release((value) w), release(twotow);
  247.     if (!Integral(result)) syserr(MESS(704, "app_floor: result not integral"));
  248.     return result;
  249. #else !EXT_RANGE
  250.     return (value) mk_int(floor(Frac(u)));
  251. #endif !EXT_RANGE
  252. }
  253.